//__gnu_debug::set
#include "pch.h"


// DEBUGGING

#define NDEBUG = 1
#define BOOST_DISABLE_ASSERTS
//#define DEBUG_STL = 1


#include <cstdint>










//
//  Standard headers
// 
#include <cassert>

#include <iostream>
using std::cerr;
using std::cout;
using std::ofstream;
using std::istream;
#include <utility>
using std::make_pair;

#include <set>
using std::set;
#include <vector>
using std::vector;
#include <string>
using std::string;
#include <numeric>
#include <algorithm>
using std::fill_n;
#include <stdexcept>
using std::runtime_error;
#include <unordered_map>
using std::unordered_map;
#include <unordered_set>
using std::unordered_set;
#include <tuple>
using std::get;
using std::tuple;
using std::tie;
using std::pair;


//
//  My common headers
// 
//#include <unique_file_name_with_replacement.h>
#include <progress_indicator.h>
#include <tuple_hash.h>
//#include <commify.h>
//#include <portable_timer.h>
//#include <intervals_per_contig.h>
//	using intervals_per_contig::t_intervals_per_contig;
//	using intervals_per_contig::t_interval_search_cache;
//	using intervals_per_contig::t_intervals;
//	using intervals_per_contig::ci_map_str_intervals;
//
//#include <open_stream_name_with_substitutions.h>
//
//
//#include <get_col_names_from_header.h>
//using parsing::parse_error;
//using parsing::get_col_names_from_header;





//
//  Boost headers
// 
//#include <boost/algorithm/string/replace.hpp>
//	using boost::replace_all_copy;
#include <boost/algorithm/string.hpp>
//    using boost::join;
using boost::algorithm::replace_all;
#include <boost/range/iterator_range_core.hpp>

///#include <boost/lexical_cast.hpp>
//	using boost::lexical_cast;
//#include <boost/regex.hpp>
//	using boost::regex;
//#include <boost/functional/hash.hpp>
//	using boost::hash_combine;
//#include <boost/iterator/counting_iterator.hpp>
//	using boost::make_counting_iterator;



//
//	Boost Spirit
// 
#include <boost/spirit/include/qi.hpp>
	namespace qi = boost::spirit::qi;
//#include <boost/spirit/include/support_multi_pass.hpp>
//#include <boost/spirit/include/karma.hpp>
//	namespace karma = boost::spirit::karma;
//	using boost::spirit::karma::uint_;
//	using boost::spirit::karma::lit;
//	using boost::spirit::karma::generate;
//#include <boost/fusion/include/adapt_struct.hpp>
//#include <boost/fusion/include/io.hpp>
//#include <boost/spirit/include/phoenix_core.hpp>
//namespace phoenix = boost::phoenix;
//#include "qi_parse_error_handler.h"
//
#include <boost/filesystem.hpp>
	using boost::filesystem::path;
	using boost::filesystem::file_size;
	using boost::filesystem::is_regular_file;
	using boost::filesystem::basename;
//
//#include <boost/array.hpp>
//	using boost::array;
//#include <boost/foreach.hpp>
//
//#include <boost/range/algorithm_ext/push_back.hpp>
//    using boost::push_back;
#include <boost/range/numeric.hpp>
	using boost::accumulate;
#include <boost/range/algorithm.hpp>
//	using boost::count;
//	using boost::sort;
//	using boost::reverse;
//	using boost::for_each;
//#include <boost/range/numeric.hpp>
//	using boost::accumulate;
//#include <boost/range/adaptors.hpp>
//
//
//#include <boost/iostreams/device/file.hpp>
//#include <boost/iostreams/stream.hpp>
#include <boost/iostreams/device/mapped_file.hpp>
	namespace io = boost::iostreams;
#include <boost/assert.hpp> 


#include <murmurhash/MurmurHash3.h>

#include <boost/format.hpp>
using boost::format;


#include "boost/multi_array.hpp"
using boost::multi_array;
using boost::const_multi_array_ref;
using boost::extents;

#ifdef DEBUG_STL
	#include <debug/set>
	#include <debug/vector>
using __gnu_debug::set;
using __gnu_debug::vector;
#else
using std::set;
using std::vector;
#endif



#include "read_options.h"


void do_main(prog_options& opt);

//________________________________________________________________________________________
// 
//  	Main
//			boilerplate to catch errors
// 
//________________________________________________________________________________________
int print_error_message(string error_string)
{
	replace_all(error_string, "\n", "\n\t");
	std::cerr
	<< "\n"<< string(60, '=') << "\n"
	<< "\nError:\n\t" << error_string 
	<< "\n\n" << string(60, '=') << "\n\n";
	return -1;
}
int main (int argc, char *argv[])
{
	try
	{
		// 
		//  Read options
		//  	return if user chose help or version
		// 
		prog_options opt;
		if (!read_options(opt, argc, argv))
			return 0;

		do_main(opt);

	}
	catch (boost::program_options::validation_error& e)
	{
		return print_error_message(e.what()  + string(" of ") + e.get_option_name());
	}
	catch (boost::program_options::multiple_occurrences& e)
	{
		return print_error_message(e.what()  + string(" of ") + e.get_option_name());
	}
	catch (const std::exception& e)
	{
		return print_error_message(e.what());
	}
	catch (...)
	{
		return print_error_message("Unexpected exception of unknown type.");
	}
	return 0;
}



void throw_logic_error(format& f)
{
	throw std::logic_error(f.str());
}



#include <boost/tuple/tuple_comparison.hpp>

typedef tuple<uint32_t, uint32_t, uint32_t, uint32_t> uint128_t;

struct hash_uint128_t
{
	size_t operator()(uint128_t const & t) const{
        return get<0>(t) + get<1>(t) + get<2>(t) + get<3>(t);
    }
};


//	sequence	type1	type2	
enum eFRAG_TYPE{eU, eF, eR};
enum eSAM_FLAGs 
{
	eMULTIPLE_FRAGS = 0x1	, // > 1 fragment
	eFIRST 			= 0x40	, // the first fragment in the template
	eLAST 			= 0x80	, // the last fragment in the template
	eALLALIGNED		= 0x2	, // all aligned
	eONLY1ALIGNED	= 0x8	, // only this fragment aligned
	eREVERSE		= 0x10	, // Reverse
	eNEXTREVERSE 	= 0x20	, // Next reverse (don't trust...)
	eUNMAPPED		= 0x4	, // unmapped: all bets off
	e2NDARY		 	= 0x100	, // secondary alignment (don't trust this alignment?)
	eQC_FAILED		= 0x200	, // not passing quality controls (don't trust this alignment?)
	eDUPLICATE		= 0x400	, // PCR or optical duplicate (don't trust this alignment?)
	eDODGY			= 0x4 | 0x100 | 0x200 | 0x400
};


struct t_frag_pair
{
	t_frag_pair(char cnt_Ns_1_, eFRAG_TYPE type_1_ = eU, char MQ_1_ = -1, int pos_1_ = -1, char cnt_Ms_1_ = -1):
					diff(0),
					type_1(type_1),
					type_2(eU),   
					MQ_1(MQ_1_),     
					MQ_2(-1),     
					cnt_Ns_1(cnt_Ns_1_),  
					cnt_Ns_2(-1),  
					cnt_Ms_1(cnt_Ns_1_),  
					cnt_Ms_2(-1),  
					pos_1(pos_1_),     
					pos_2(-1){}
	int				diff;
	eFRAG_TYPE 		type_1;
	eFRAG_TYPE 		type_2;
	char			MQ_1;
	char			MQ_2;
	char			cnt_Ns_1;
	char			cnt_Ns_2;
	char			cnt_Ms_1;
	char			cnt_Ms_2;
	int 			pos_1;
	int 			pos_2;
};


typedef std::unordered_map<uint128_t, t_frag_pair, hash_uint128_t> map_llfp;
typedef std::unordered_map<uint128_t, unsigned, hash_uint128_t> map_llu;

void initial_calls_survey(	prog_options& 					opt, 
							const string&					file_name,
							map_llfp&						reads_per_template,
							map_llu&						multiples_hits)
{
	//BOOST_FOREACH(const string& s, contig_names) 
	opt.both() << format("Parsing %s\n") % file_name;
	
	
	// initialise frag_pairs_to_names for output	
	unordered_map<int, char> frag_type_to_names {{eU, 'U'}, {eF, 'F'}, {eR, 'R'}};
	//unordered_map<tuple<int, int>, string, tuple_hash<tuple<int, int>>> frag_pairs_to_names;
	unordered_map<tuple<int, int>, string> frag_pairs_to_names;
	auto frags_types = {eU, eF, eR};
	BOOST_FOREACH(auto ii, frags_types) 
		BOOST_FOREACH(auto jj, frags_types)
			frag_pairs_to_names[std::make_tuple(int(ii), int(jj))] = frag_type_to_names[ii] + frag_type_to_names[jj];

	// 
	// guess number of rows for progress dots based on file sizes
	// 
	//path input_path(file_name);
	if (!is_regular_file(file_name))
		throw std::logic_error(file_name + " does not exist or is not a regular file.");

	// 17 bytes per line for progress dots
	auto input_file_size = file_size(file_name);

	io::stream<io::mapped_file_source> input_stream;
	input_stream.open(file_name.c_str());
	if (!input_stream)
		throw std::logic_error("Can't open " + file_name);


	// 
	//	Parsing grammar 
	// 
	using qi::graph;
	using qi::uint_;
	using qi::int_;
	using qi::lit;
	using qi::char_;
	using qi::omit;
	using qi::raw;
	auto space = qi::omit[qi::space];
	auto grammar = 		raw[+graph]		> space 		//	0 QNAME
					>	uint_			> space   		//	1 FLAG  
					>	omit[+graph]	> space   		//    RNAME 	Contig
					>	uint_			> space   		//  2 POS   
					>	uint_			> space   		//  3 MAPQ  	

					// 	Matching counts from cigar 		//  4 CIGAR
					//  \*|([0-9]+[MIDNSHPX=])+
					//	Strictly speaking this field has an asterisk OR cigar string
					//		our grammar has (optional) asterisk followed by cigar string
					//		so that we always have vector<int> rather optional
					>	omit[-lit('*')] >				
						*(	omit[*(uint_ >> char_("IDNSHPX="))] >>
							uint_ >> lit('M') >>
							omit[*(uint_ >> char_("IDNSHPX="))]) > space	

					>	omit[+graph]	> space   		//    RNEXT 
					>	uint_			> space   		//  5 PNEXT	Position of next 
					>	omit[int_] 		> space   		//    TLEN  
					>	raw[+graph]		> space   		//  6 SEQ   
					>	omit[+graph]	> -space   		//    QUAL  
					> 	omit[-(+graph)];				//	  any other stuff
					;



	unsigned cnt_rows = 0;
	string line_str;

	// Two rows of dots progress dots
	auto bytes_read				= input_file_size * 0;

	try
	{
		unsigned guess_rows_length	= 250;
		t_progress_indicator progress(opt.verbose(), input_file_size, guess_rows_length, 100);
		progress.use_timer();
		progress.set_threshold_for_completion_time_estimate(1);
	
		while (getline(input_stream, line_str))
		{
			bytes_read += line_str.length();
			progress.increment(input_file_size, bytes_read, cnt_rows++);

			if (!line_str.length() || line_str[0] == '@')
				continue;



			//
			// 		Parse
			// 
			typedef boost::iterator_range<string::iterator> string_range;
			string_range		qname;				// 0
			unsigned			flag;   			// 1
			unsigned			pos;				// 2
			unsigned			mapq;   			// 3
			vector<unsigned>	cigar_matches;  	// 4
			unsigned			pnext;  			// 5
			string_range		sequence;			// 6

			auto beg = line_str.begin();
			auto end = line_str.end();
			qi::parse(beg, end, grammar, qname, flag, pos, mapq, cigar_matches, pnext, sequence);

			// count matches
			unsigned cnt_matches = accumulate(cigar_matches, 0U);
			// count Ns
			unsigned cnt_Ns = boost::count(sequence, 'N');

			
			uint128_t qname_hash;
			MurmurHash3_x86_32 (&qname.front(), qname.end() - qname.begin(), 0, &qname_hash);

			if (flag & eDODGY)
			{
				auto inserted = reads_per_template.insert({qname_hash, t_frag_pair(cnt_Ns)});
				if (!inserted.second)
				{
					if (inserted.first->second.cnt_Ns_2 >= 0)
						++multiples_hits[qname_hash];
					else
						inserted.first->second.cnt_Ns_2 = cnt_Ns;
				}
			}

			// non-dodgy: insert into first or second position unless is third time, in which case
			// ++ multiple
			else
			{
				eFRAG_TYPE frag_type = flag & 0x10 ? eR : eF;

				auto inserted = reads_per_template.insert({qname_hash, t_frag_pair(cnt_Ns, frag_type, mapq, pos, cnt_matches)});
				if (!inserted.second)
				{
					auto& pp = inserted.first->second;
					if (pp.cnt_Ns_2 >= 0)
						++multiples_hits[qname_hash];
					else
					{
						pp.cnt_Ns_2 = cnt_Ns;
						pp.type_2	= frag_type;
						pp.MQ_2		= mapq;
						pp.pos_2	= pos;
						pp.cnt_Ms_2 = cnt_matches;
					}
				}
			}





			//if (cnt_rows > 100)
			//	break;

		}
	}
	catch (qi::expectation_failure<string::iterator>& e)
	{
		throw_logic_error(format("Problems parsing line #%d from \"%s\" (%s at <%s>):\n[%s]") 
			% cnt_rows % path(file_name).filename().native() % e.what() % string(e.first, e.last) % line_str);
	}
}



void do_main(prog_options& opt)
{
	map_llfp						reads_per_template;
	map_llu							multiples_hits;
	BOOST_FOREACH(const string& input_filename, opt.input_filenames) 
	{
		initial_calls_survey(	opt, input_filename, reads_per_template, multiples_hits);
	}

	unordered_map<unsigned, unsigned> cnts;
	BOOST_FOREACH(auto& pp, multiples_hits) 
	{
		cnts[pp.second] ++;
	}

	BOOST_FOREACH(auto& pp, cnts) 
	{
		opt.both() << format("%1$d items with %2$d reads\n") % pp.second % (pp.first + 2);
	}


}



